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Final  Report 

ABSTRACT 

Several  classes  of  algorithms  for  solution  of  the  general 
nonlinear  programming  (constrained  optimization)  problem,  and  four 
specific  implementations  of  these  were  chosen  and  evaluated  with 
respect  to  expected  speed  of  computation.  A  test  problem  based  on  the 
path  generation  problem  of  terrain  following  /  terrain  avoidance 
flight  was  developed,  and  the  performance  of  the  chosen  optimization 
procedures  was  compared.  It  was  found  that  the  generalized  reduced 
gradient  method  was  faster  and  more  reliable  than  either  of  two 
augmented  Lagrangian  methods  and  a  quadratic  approximation  method. 
However,  the  solution  time  for  the  TF/TA  type  problem  was  found  to  be 
far  in  excess  of  what  would  be  required.  Several  simplifications  of 
the  problem  statement  were  attempted  in  order  to  decrease  computation 
time  without  compromising  the  integrity  of  the  solution. 
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I.  INTRODUCTION 

While  participating  as  a  research  associate  in  the  Air  Force  Summer 
Faculty  Research  Program  as  part  of  the  Avionics  Laboratory  at  Wright- 
Patterson  AFB,  the  author  became  involved  in  the  simulation  of  automatic 
terrain  following/terrain  avoidance  (TF/TA)  flight.  It  became  apparent 
that  the  "path  generation  problem",  as  currently  stated,  involves  a  huge 
computational  burden  due  to  the  large  number  of  variables  and  the  large 
number  of  state  and  control  constraints.  While  in  theory  the  problem 
can  be  solved  using  known  techniques,  the  computation  time  necessary  to 
solve  it  is  far  beyond  the  time  in  which  the  solution  is  needed.  While 
the  development  of  very  high  speed  integrated  circuit  technology  may  lessen 
this  gap,  it  appears  that  new  computational  procedures  which  decrease  the 
number  of  calculations  for  constrained  optimization  will  be  necessary. 

The  overall  thrust  of  the  research  project  was  to  study  and  evaluate 
current  nonlinear  optimization  algorithms  as  to  their  relative  success  in 
dealing  with  large  problems  and  to  suggest  variations  in  the  methods  which 
could  lead  to  improved  performance.  Computation  time  was  the  key  variable 
of  interest  in  the  study.  After  evaluations  were  completed  on  a  test  set 
of  general  large  problems  and  in  particular  on  a  TF/TA  test  problem,  atten¬ 
tion  was  paid  to  assessing  several  simplifications  of  the  TF/TA  path 
generation  problem  since  that  was  one  of  the  driving  forces  for  the  invest¬ 
igation.  For  example,  in  the  TF/TA  problem  it  is  important  to  have  "optimal" 
points  for  the  beginning  of  the  path,  while  it  may  only  be  necessary  to 
obtain  approximate  values  for  the  rest  of  the  segment. 

A  first  research  objective  was  to  evaluate  candidate  algorithms  from 
within  the  general  classes: 

a)  generalized  reduced  gradient  techniques 
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b)  feasible  directions  techniques 

c)  augmented  Laqrangian  techniques  (multiplier  methods) 
and  d)  quadratic  approximation  methods 

as  to  their  success  in  solving  large  (in  terms  of  both  number  of  variables 
and  number  of  constraints)  nonlinear  optimization  problems.  Particular 
implementations  of  these  algorithms  differ  in  the  mechanisms  used  for 
carrying  out  line  searches,  unconstrained  minimizations,  nonlinear  equation 
solving,  etc..  Efficient  computation  through  the  use  of  such  procedures 
as  quasi-Newton  updating,  quadratic  approximation,  and  parallel  processing 
was  addressed.  Several  specific  implementations  (computer  code)  for  the 
various  algorithms  were  obtained. 

We  then  classified  the  path  generation  problem  as  to  the  number  of 
variables,  number  of  linear  constraints,  number  and  forms  of  the  nonlinear 
constraints,  bounds  on  the  variables,  etc.,  and  selected  reasonable  values 
for  the  constraints  in  a  test  problem.  The  mathematical  programming  test 
problem  was  scaled  in  order  to  improve  its  numerical  properties. 

Each  specific  implementation  of  a  computational  method  had  its  own 
input  format  and  exact  problem  formulation,  and  so  the  test  problem  was 
restated  in  forms  suitable  for  solution  by  each  of  these. 

A  model  for  terrain  data  was  developed  and  this  provided  realistic 
terrain  data  to  the  test  problem.  Several  stylized  terrains  were  also 
programmed  in  order  to  test  the  performance  of  the  algorithms. 

The  major  project  effort  was  the  comparison  and  evaluation  of  several 
of  the  "better"  methods  for  constrained  nonlinear  optimization,  with  the 
path  generation  problem  used  as  a  test. 

Finally,  the  effects  that  some  alternative  statements  of  the  optim¬ 
ization  problem  have  on  solution  time  and  the  accomplishment  of  the  TF/TA 
objective  were  addressed. 


In  Section  II  of  this  report  the  background  and  mathematical  state¬ 
ment  of  the  terrain  following/  terrain  avoidance  path  generation  problem 
is  given.  In  Section  III  a  model  for  terrain  data  is  developed.  Inter¬ 
polation  procedures  for  this  data  are  also  discussed. 

A  survey  of  numerical  methods  for  general  nonlinear  programming 
problems  appears  in  Section  IV,  while  several  particular  implementations 
of  the  algorithms  are  described  in  Section  V.  In  Section  VI  we  state 
our  test  problem,  its  parameters,  and  computational  experiments  that 
were  carried  out.  The  results  of  these  experiments  are  stated  and 
discussed.  Also,  several  modifications  of  the  original  problem  are 
performed  and  these  are  evaluated. 

Section  VII  restates  project  results  and  conclusions.  This  is 
followed  by  a  Bibliography  which  includes  a  survey  of  relevant  literature 
and  all  references  that  are  cited  in  the  text  of  the  report. 

The  duration  of  the  project  was  three  months  and  it  required  the 
one-half  time  dedication  of  the  principal  investigator  during  this  period. 
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II.  THE  TERRAIN  FOLLOWING  /  TERRAIN  AVOIDANCE  PROBLEM 

During  a  ten  week  period  in  the  summer  of  1981,  the  author  developed 
a  simulation  tool  for  the  evaluation  of  terrain  following/terrain  avoidance 
flight  while  working  with  the  System  Concepts  group  of  the  Mission  Avionics 
division  at  the  Avionics  Laboratory  of  the  Wright  Aeronautical  Laboratories 
at  Wright-Patterson  AFB,  Ohio.  Current  approaches  to  terrain  following  and 
terrain  avoidance'  '  were  studied  and  models  for  the  components  of  a  TF/TA 
system  were  developed.  A  final  report  was  written  which  documented  the 
results  obtained  during  this  period  [ref.  1].  While  working  on  the  simulation 
package  it  became  apparent  to  the  author  that  proposed  methods  of  calculating 
desired  flight  paths  would  not  be  able  to  perform  the  calculations  in  near 
real  time,  and  that  improvement  will  be  necessary  before  the  process  can 
be  implemented.  A  description  of  the  TF/TA  problem,  its  importance,  and 
methods  of  solution  will  be  given  at  this  point. 

Sustained  low  level  flight  is  one  key  to  successful  penetration  into 
a  closely  monitored  region.  By  flying  close  to  the  ground  the  probability 
of  detection  by  ground  based  or  air  based  radar  systems  is  substantially 
reduced  due  to: 

1'  direct  masking  of  the  vehicle  by  the  terrain  between  it  and 
the  source  of  the  illumunation  (Fig.  la) 

2)  indirect  masking  due  to  the  inability  of  the  radar  system  to 
distinguish  the  vehicle  from  ground  clutter.  (E.  G. -vehicle  and 
terrain  are  within  the  same  range  or  angle  resolution  cell).  (Fig.  lb,  lc) 


(*)  For  the  purposes  of  this  proposal  the  term  "terrain  following"  refers 
to  vertical  maneuvers  (go  o •  wh  a  "terrain  avoidance"  refers  to  horizontal 
and  vertical  maneuvering  (go  around). 


(S  aircraft 


Figure  la 
Direct  Masking 


Figure  lb 
Indirect  masking 


Figure  1c 
Indirect  Masking 


Even  if  detection  is  accomplished,  low  level  flight  decreases  the  probability 
of  a  successful  track  being  initiated. 

On  the  other  hand,  the  probability  of  collision  with  the  ground  cr 
with  obstacles  (houses,  towers,  transmission  lines,  etc.)  is  substantially 
increased  as  lower  clearance  heights  are  attempted.  In  order  to  avoid 
collision  it  is  necessary  that  an  accurate  description  of  present  and 
upcoming  terrain  features  and  obstacles  be  available  to  the  flight  path 
generator  and  flight  control  system  so  that  a  safe  path  can  be  flown. 

An  approach  that  has  been  taken  for  automatic  terrain  following  is 
depicted  below  (Fig.  2). 


A  forward  looking  radar  scans  in  elevation  from  an  angle  ^  to  . 

The  return  signal  provides  range  and  angle  of  elevation  information  of 
the  terrain  with  respect  to  the  aircraft.  A  flight  path  is  chosen  that 
will  allow  clearance  of  all  terrain  points  by  some  predetermined  height. 
The  development  of  the  algorithm  which  specifies  flight  commands  based 
on  such  information  is  described  in  the  series  of  ADLAT  studies  carried 


out  by  Calspan  Inc.  [ref.  2,3].  A  radar  based  automatic  terrain  following 
system  is  operational  in  the  Air  Force  F-111B  aircraft.  This  system  is 
designed  for  a  minimum  clearance  height  on  the  order  of  a  few  hundred  feet. 

Unfortunately  this  type  of  system  exhibits  the  property  known  as 
"ballooning"  as  depicted  in  Figure  3. 


The  reason  for  ballooning  is  that  the  system  is  unable  to  see  parts  of  the 
terrain  that  are  masked  from  its  line  of  sight  (Figure  4). 


For  terrain  following  and  terrain  avoidance  a  scan  in  both  azimuth  and 
elevation  would  be  used,  an  algorithm  developed  to  specify  flight  commands, 
and  then  these  would  be  implemented  by  the  flight  control  system.  Terrain 
masking  would  remain  a  problem  and  would  severely  limit  side  to  side 
maneuvers  since  ballooning  in  the  horizontal  plane  could  prove  disastrous. 

Use  of  the  Digital  Land  Mass  Simulation  (DLMS)  data  base  from  the 
Defense  Mapping  Agency  has  recently  been  considered  as  a  potential  source 
of  information  to  a  terrain  following/terrain  avoidance  system.  The  DLMS 
Level  II  data  base  breaks  the  surface  of  the  earth  into  approximately 
100  foot  square  grids  and  provides  smoothed  altitude  values  for  each  of 
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these  grid  areas.  Incorporation  of  this  data  into  a  terrain  following 
or  terrain  avoidance  system  as  a  stored  on-board  map  allows  greater 
look-ahead  and  fills  in  areas  that  are  masked  from  the  sensors.  The  use 
of  map  data  for  the  TF/TA  application  involves  two  difficulties.  First, 
an  accurate  fix  of  the  aircraft  position  is  required  so  that  the  correct 
map  portion  is  called  down  from  memory.  Second,  the  DLMS  data  may  not 
include  obstacle  information. 

An  accurate  fix  of  aircraft  position  might  be  obtained  by  periodically 
updating  the  inertial  navigation  system  (INS)  by  either  terrain  contour 
correlation  (comparison  of  altimeter  data  to  a  stored  map)  or  from  the  ‘ 

Global  Positioning  System  (GPS).  Once  a  position  is  calculated,  a  section 
of  the  DLMS  data  would  be  called  down  and  used  to  augment  sensor  data  in 
order  to  provide  a  "best  guess"  for  the  terrain  ahead.  ; 

Work  is  presently  underway  to  address  the  three-dimensional  terrain 
avoidance  problem,  with  the  goal  of  increased  survivability  and  lower  minimum 
clearance  height  [ref.  4,5],  Key  to  the  success  of  this  work  will  be  the 
necessity  to  accurately  characterize  what  is  ahead  and  to  choose  a  best 

path  based  on  this  knowledge.  I 

A  depiction  of  the  TF/TA/OA  system  is  shown  in  Figure  5.  f 


(choosaa  •  path  for  the  naxt 
t j  •  S  t  jaconda) 

Figure  5 


The  path  generator  is  the  heart  of  the  system.  It  searches  the 
augmented  map  and  prescribes  a  path  that  attains  low  altitude  flight 
within  the  constraints  of  aircraft  performance  and  crew  comfort.  Lateral 
deviation  from  a  prespecified  ground  track  may  also  be  weighted  in  the 
optimization.  The  flight  control  system  then  provides  control  signals  to 
the  actuators  for  the  aerodynamic  surfaces  and  attempts  to  follow  this 
optimal  path.  It  is  anticipated  that  the  desired  flight  path  would  be 
computed  for  several  segments  ahead  and  recomputed  as  each  segment  is 
flown  and  new  sensor  data  is  obtained.  The  update  time  t  would  be 
determined  by  such  factors  as  terrain  roughness,  speed  of  the  aircraft, 
and  ranges  of  the  sensors.  The  calculation  of  the  desired  path  must  be 
performed  within  t  seconds.  It  should  also  be  noted  that  since  only 
ty  seconds  of  an  N-tu  second  flight  path  will  be  flown,  it  is  only 
necessary  that  the  beginning  of  the  calculated  best  path  be  optimal, 
provided  that  no  constraints  are  violated  on  the  rest  of  the  path.  This 
is  for  protection  against  the  case  where  optimization  in  the  short  run 
leads  to  the  necessity  to  perform  impossible  maneuvers  to  avoid  a  crash. 

The  time  N't^  must  be  long  enough  so  that  a  maneuver  around  a  worst  case 
terrain  point  will  be  possible. 

A  diagram  depicting  the  sequence  of  events  is  shown  in  Figure  6.  It 
is  expected  that  an  update  time  of  1-2  seconds  and  N*tu  =  t^  of  5-10  seconds 
will  be  necessary  for  adequate  performance  of  the  system.  In  essence  this 


Figure  6" 
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means  that  a  5  second  long  flight  path  must  be  generated  by  the  path 
generator  each  second. 

As  described  above,  the  TF/TA  optimization  procedure  can  be  seen  as 
that  of  minimizing  a  function  of  altitude  and  lateral  deviation  (cost 
function)  under  a  set  of  constraints  given  by 

1)  aircraft  state  equations 

2)  aircraft  performance  limits  and  crew  comfort  restrictions 

(state  and  control  constraints) 

3)  terrain  restrictions  (trajectory  constraints). 

For  a  simple  quadratic  cost  functional  and  a  point  mass  model  for  the 
aircraft  the  optimization  problem  has  been  stated  in  [ref.  5]  as: 

For  each  path  segment  of  length  N*t  =  t^  seconds 

min 

P(t),  nz(t)  0 

subject  to  the  state  equations  (point  model  for  aircraft) 


ff  { pz ^  +  •  [ Py ( t )  -  yd(t)]  -dt 


Px  =  V  cos  t cos  ^ 
p y  =  V  cos'/sin 


three  dimensional 
position 


Pz  =  V  sin  /  J 

-  -(nzsin<^)/ (V  cos  ^  )  heading  angle 

* 

7  =  (nz  cosf  -  g  cos/)  /  V  flight  path  angle 
ft  =  p(t)  roll  rate 


V  =  -g  sin  Y 


roll  rate 


velocity 


and  at  each  trajectory  point  the  constraints 

Pz(t)  ^  terrain  height  plus  safety  factor 


7.  -  V(t)  -  * 

min  ma 

fmin  “fm 


flight  path  angle 


bank  angle 
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«i  *  * 

<p  .  i  <$U)  - 

T  mm  ~  v  '  7  max 


roll  acceleration 


/t>.  ^/n(t)  -  load  factor  (determined  by 

zmin  *  zmax  allowable  ride  level) 

/£  AA(t)-  &  pitch  jerk 

min”  e  max 

The  control  variables  for  the  problem  are  n2(t)  (load  factor)  and  p(t) 
(roll  rate). 

In  generality  the  problem  is  of  the  form: 

min  fQ(x  ,  u  ,  t  ) 
u^ 

subject  to  a  £  x  £  b  state  constraints 

c  4=  u  is  d  control  constraints. 


As  it  is  stated  above,  the  problem  is  a  continuous  time  nonlinear 
optimal  control  problem  with  state  and  control  constraints.  The  application 
of  the  variational  approach,  the  Pontryagin  Maximum  principle,  or  continuous 
time  dynamic  programming  leads  to  computationally  infeasible  methods  of 
solution  [ref.  6]. 

Discretization  in  a  time  step  dt  transforms  the  original  infinite 
dimensional  problem  (find  the  functions  nz(t)  and  p(t))j  into  a  finite 
dimensional  one  (find  sequences  n..(0), - »n  (N-l),  p(0), _ ,p(N-l ) ) . 

The  problem  becomes:  N 

Problem  PI:  min  ^/pz(i)  +  W  [py(i )  -  yd(i )]  f 

p(Q)  - - >p(N-l )  1=1 

n2(0), - ,nz(N-l) 

subject  to  the  state  equations 

x(i+l)  =  f  (x(i),  nz(i),  p(i)  )  for  all  i=l,..,N-l 

and  a  set  of  constraints  which  most  hold  for  each  i=l,..,N-l. 
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The  complexity  of  this  optimization  problem  clearly  depends  on  the 
size  of  the  time  step  At  and  the  final  time  t^.  A  small  discretization 
interval  leads  to  a  good  approximation  to  the  continuous  time  system  and 
avoids  violation  of  constraints  due  to  too  coarse  a  grid,  but  it  creates 
a  larger  optimization  problem. 

There  are  several  different  ways  of  looking  at  the  problem  PI.  First, 
it  can  be  viewed  as  a  discrete-time  optimal  control  problem  subject  to  the 
state  and  control  constraints.  Invocation  of  the  discrete  maximum  principle 
[ref.  6,7]  recasts  the  solution  of  the  optimization  as  the  solution  of  a 
two  point  boundary  value  problem.  Because  of  the  state  constraints  and 
since  the  state  equations  are  nonlinear,  the  solution  of  the  resulting 
boundary  value  problem  proves  difficult  to  obtain. 

A  second  approach  views  the  optimization  problem  as  an  N  stage 
sequential  decision  process.  Adopting  this  point  of  view  the  dynamic 
programming  algorithm  [ref.  6,  8,  9,  10,  11,  12]  can  be  used  to  construct 
the  optimal  control  sequences  and  optimal  trajectory.  Under  this  approach 
the  large  number  of  constraints  actually  lessen  the  amount  of  computation 
that  would  be  necessary  for  the  unconstrained  case. 

The  basis  for  dynamic  programming  is  the  "Bellman  Principle  of 
Optimality"  which  can  be  stated: 

"  an  optimal  policy  has  the  property  that,  whatever  the  initial 
state  and  initial  decision  are,  all  remaining  decisions  must  constitute 
an  optimal  policy  with  regard  to  the  state  resulting  from  the  first 
decision". 

The  principle  can  be  applied  to  multi-stage  decision  problems  in  that 
it  can  be  seen  that  the  minimum  cost  from  state  x  at  a  stage  k  is  found 
by  minimizing  the  sum  of  the  current  single  stage  cost  plus  the  minimum 
cost  of  going  to  the  end  of  the  process  from  the  resulting  next  state. 
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Unfortunately  for  all  but  the  simplest  problems  the  dynamic  programming 
algorithm  is  impractical  due  to  what  has  become  to  be  known  as  "the  curse 
of  dimensionality".  In  order  to  evaluate  minimum  costs  at  each  stage  it 
is  necessary  to  discretize  the  state  and  control  variables  so  that  they 
may  take  on  only  a  limited  number  of  values.  For  example,  a  roll  angle 
might  be  discretized  to  the  set  £-30°,  -20°,  -10°,  0°,  10°,  20°,  30 °J  .  If 
for  example,  there  are  7  state  variables,  each  discretized  into  only  10 
values  the  dynamic  programming  solution  requires  a  high  speed  memory  of 
107  words,  which  is  beyond  the  capacity  of  most  large  systems. 

Computation  time  is  a  function  of  the  number  of  stages,  and  for  a  20 
stage  problem  there  would  be  on  the  order  of  20  x  107  calculations  required, 
and  at  a  rate  of  1  million  calculations  per  second  this  would  take  about 
200  seconds.  Another  disadvantage  of  the  dynamic  programming  approach  is 
that  no  part  of  the  optimal  trajectory  is  found  until  the  complete  trajectory 
is  calculated.  The  previously  discussed  factors  make  it  clear  that  this 
procedure  is  not  suitable  for  the  TF/TA  problem  as  stated. 

A  third  way  of  viewing  Problem  PI  is  as  a  general  constrained  nonlinear 
optimization  problem.  Under  this  point  of  view  the  problem  is: 

min  f(  ^  ) 
z 

subject  to  g.j(z)^  0  i=l,..,m  inequality  const. 

h  • (zj  -  0  j=l,...,q  equality  const. 

J 

where  the  vector  z  contains  all  control  values  and  state  variable 
values  for  all  time  steps.  That  is: 


z  =  [nz(0), . nz(N-l),p(0), - ,p(N-l ) ,px(0) , . . . ,PX(N-1 ) ,py(0) , . 

Py ( N - 1 ) . V ( 0 ) , . . . ,  V ( N-l )  ]t. 

The  functions  h.(z)  2  0  include  the  state  equation  relations  at 

J 
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each  1=1, . ,N,  while  the  g^(z)%0  include  all  inequality  constraints 

on  state  and  control  at  each  stage  i = 1 . ,N. 

It  is  seen  that  this  is  an  optimization  over  9-N  variables  subject 
to  a  large  set  of  equality  and  inequality  constraints.  Although  this 
problem  is  formidable,  several  methods  have  been  developed  for  the  general 
nonlinear  programming  problem.  An  outline  of  some  of  the  more  recent  and 
successful  methods  appears  in  Section  IV. 


14 


III.  TERRAIN  MODEL  AND  INTERPOLATION 


One  of  the  TF/TA  constraints  was  that  the  altitude  of  the  airplane 
be  greater  than  a  fixed  distance  above  the  terrain.  This  constraint  is 
of  the  form 

Z  -  h(X,Y)»  Znjn 

where  Z  is  the  altitude,  and  h(X,Y)  is  the  terrain  height  at  position 
(X,Y).  In  order  to  develop  a  realistic  test  problem,  and  since  it  is 
anticipated  that  some  form  of  stored  and  updated  map  would  be  used  we 
decided  to  make  a  discrete  terrain  model.  The  model  has  adjustable  rough¬ 
ness  parameters  in  order  to  test  the  optimization  procedure  for  different 
terrain  types.  In  this  section  we  first  describe  the  terrain  model  and 
then  proceed  to  discuss  methods  for  approximating  data  at  non-grid  points. 

A.  Terrain  Model 

Statistical  properties  of  mean,  trend  removed,  actual  terrain  have 
been  shown  to  be  fairly  well  characterized  as  Gaussian  distributed,  with 
exponentially  decaying  correlation  among  samples.  The  exponential  correlation 
function  is  of  the  form 

R(r)  -  e|z(x)Z(x+7’)J  =  A  expC-T* /Tauc  ) 

where  the  factor  Tauc  is  called  the  decorrelation  distance.  Decorrelation 
distances  for  typical  terrain  samples  were  found  to  vary  from  2500ft.  (fairly 
rough  terrain)  to  30,4)00  ft.  (fairly  smooth  terrain).  An  illustration  of  the 
effect  of  Tauc  is  given  in  Figure  7.  Both  cases  would  have  the  same 
mean  value  and  standard  deviation. 


Short  decorrelation 
distance 


B 


Long  decorrelation 
distance 
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For  the  purposes  of  the  simulation  package  it  was  decided  to  produce 
synthetic  terrain  as  "real  world".  Therefore  it  was  necessary  to  develop 
a  procedure  to  generate  points  which  possessed  the  desired  statistical 
properties.  These  properties  were: 

1)  Terrain  characterized  as  a  two-dimensional  stochastic 
process  Z(x,y). 

2)  Heights  of  the  terrain  at  each  (x,y)  point  are  Gaussian 
distributed  with  mean  value  Terrmn,  and  standard  deviation 
Sdterr. 

3)  Terrain  heights  correlated  in  both  dimensions  with  an 
exponential  correlation  function.  The  decorrelation  distance 
was  Tauc. 

Time  series  analysis  provides  a  mechanism  for  generating  such  a 
process.  For  the  one-dimensional  case  it  is  known  (Re“f.  37  )  that  a 
first-order  autoregressive  process  of  the  form 


zt  =  a  Zt_i  +  $  +  Wt-1 

(where  Wt  is  a  zero  mean,  Gaussian  white  noise  process  of  standard  deviation 
given  by  O' ^  and  £  is  a  constant  trend  factor)  produces  the  process,  Zt. 


This  process  is  also  Gaussian,  with  a  mean  value 

1  - 

I  ~<X 


and  standard  deviation. 


I  -  «- 


Evaluation  of  the  autocorrelation  function  for  this  process  shows  that 

•f 

R  ( 7 ')  -  A  a  . 

z 

For  the  case  at  hand,  a  generalization  to  a  two-dimensional  process 
must  be  made.  (Pi 9-  ®) 


! 
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Consider  a  general  two-dimensional  first  order  discrete  autoregressive 
process  of  the  form 

Z(x,y)  =  a1  Z(x-l.y)  +  a2  Z(x,y-1)  +  a 3  Z(x-l,y-l)  +  <T  +  w(x,y)  (1) 

where  w(x,y)  is  a  Gaussian  white  noise  process  with  standard  deviation 

given  by  and  zero  mean, 
w 

It  is  desired  to  find  if  constants  a  ,  a„,  a.,*^  ,  and  exist  so 

1  J  w 

that  the  process  Z(x,y)  possesses  a  Gaussian  distribution  with  specified 
mean,  standard  deviation,  and  exponential  correlation  function  of  given 
decorrelation  distance.  By  analogy  to  the  one-dimensional  case,  and  since 
the  correlation  in  the  x  and  y  directions  should  be  the  same,  we  choose 
a-^  =  a2  =  a.  That  is 

Z(x,y)  =  a  Z(x-l,y)  +  a  Z(x,y-1)  +  b  Z(x-l,y-l)  +  +  w(x,y)  (2) 


The  autocorrelation  function  for  a  two-dimensional  discrete  process  is 
7f±.  =  E^Z(x,y)-Z(x+i,y+j)j  (3) 

and  we  desire  that  .  .  be  of  the  form 

ij 

=  (Sdterr)^.  exp(-(i+j) /Tauc)  .  (4) 

Evaluation  of  the  functions  ^00,  Pol ,  P'10,  and  ,  along  with 
algebraic  mainipulations  and  the  enforcement  that  ^io  =  /  01  and  ^  ^ 


generates  the  set  of  relations 

2!  -  (1  -  b  -  2g2  ){rl 

00  (1  -  b  -  4a^  -  4  a^  b  -  b^  +  b^) 


(5) 


By  setting  b 


/ 10 


2 

*  -a 


we  arrive  at  the  desired  relations 
2 

<T~  w 


1 1  2.2 

(1  -  a  ) 


(6) 

(7) 


(8) 
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1  LG 


*01  =  a^' 


a  0  qo 


and 


11 

Z  = 


2  )f 
a  *00 

s 

(1  -  a)2 


* 


ij 


=  a 


i+j  } 


By  setting 

a 

and  „ 

<T2 

w 


exp(-l/Tauc) 

2 

Terrmn • (1-a) 
(Sdterr)2(l-a2)2 


(9) 

(10) 
01) 


we  obtain  a  two-dimensional  stochastic  process  such  that  the  mean  is  Terrmn, 
the  standard  deviation  is  Sdterr,  and  the  process  is  exponentially 
correlated  with  decorrelation  factor,  Tauc. 

Tills  process  needs  initial  conditions.  These  are  provided  as  illustrated 
in  Figure  9. 
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Number  of 
y-points,  MPTS 


Figure  9, 


with  Z(l,l) 

ZCl.j) 

Z(i,l) 

Z(i,j) 


a  Gaussian  random  variable  with  standard 
deviation  =*  Sdterr,  and  mean  *  Terrmn. 

a  Z(l, J-l)  +  /,  +  Wj^  first  Vow 

a  Z(i-1 , 1)  +  <T,  +  w^  first  column 

a  Z(i-1, j)  +  a  Z(i, j-l)  -  a2  Z(i-l.J-l)  +  <£  + 
^  i  -  2,3,...,  NPTS 
J  -  2,3,...,  MPTS 
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Due  to  the  large  number  of  Gaussian  random  variables  that  will  be 
generated  throughout  the  simulation  a  fast  procedure  to  generate  such 
numbers  was  needed.  The  procedure  chosen  was  to  generate  a  uniformly 
distributed  random  variable  in  (0,1)  and  then  to  refer  to  a  tabulation  of 
the  inverse  normal  distribution  function  to  obtain  a  Gaussian  zero  mean, 
standard  deviation  =  1,  random  variable.  This  value  was  then  multiplied 
by  the  desired  <Tand  then  added  to  the  desired  mean  value.  This  type 
of  procedure  is  several  times  faster  than  the  more  commonly  used  procedure 
of  generating  and  adding  a  series  of  uniformly  distributed  numbers  and  then 
relying  on  the  laws  of  large  numbers. 

A  computer  program  was  written  to  generate  artificial  terrain  according 
to  the  model  developed  in  equations  (2)  through  (11).  The  input  parameters 
of  this  program  are 

Sdterr  =  the  desired  terrain  standard  deviation  in  altitude 
Terrmn  =  the  desired  terrain  mean  value  of  altitude 
Tauc  *  the  decorrelation  distance  of  the  terrain  sample 
(short  Tauc  =  rough,  long  Tauc  =  smooth) 

NPTS  and  MPTS  =  the  number  of  x  and  y  values  to  be  generated. 

The  next  several  figures  (T.l  through  T.5)  show  typical  terrain  maps 
generated  by  the  terrain  generator  program.  These  plots  were  made  using 
the  DISSPLA  package  of  plotting  routines.  It  can  be  seen  from  the  figures 
that  a  wide  variety  of  terrain  types  can  be  generated  through  variation  of 
the  input  parameters  to  the  procedure. 


ilw i 


u-WTMr.a-. 


21 


B.  Interpolation  and  Smoothing  of  the  Terrain  Data 


As  described,  the  terrain  model  generates  terrain  heights 
h(xi  •  yj>  =  hij 

for  a  discrete  grid  of  (x,y)  values.  However,  in  the  course  of  the  optim¬ 
ization  of  the  flight  path  values  of  terrain  height  at  positions  other 
than  grid  points  are  needed.  One  method  for  estimating  unknown  data 
is  to  construct  a  function  or  sequence  of  functions  which  comes  close  to 
describing  the  known  data,  and  to  use  these  same  functions  evaluated  at 


non-grid  points  as  estimates  of  missing  data.  Interpolation,  least 
squares  approximation,  and  spline  analysis  are  examples  of  this  approach. 

For  the  problem  at  hand  two  factors  must  be  considered.  First,  the 
data  is  two-dimensional  and  made  up  of  many  grid  points.  Calculating  one 
function  to  describe  all  terrain  heights  is  infeasible,  and  also  would 
lead  to  highly  oscillatory  interpolating  functions.  A  piecewise  approxi¬ 
mation  covering  a  neighborhood  of  points  is  therefore  chosen.  Second, 
since  computation  time  is  the  focus  of  the  project,  the  amount  of  compu¬ 
tation  necessary  to  arrive  at  an  approximation  should  be  as  small  as 
possible. 

Perhaps  the  simplest  procedure  for  piecewise  approximation  of  two 
dimensional  data  is  to  first  locate  the  four  closest  grid  points  to  the 
given  (x,y)  point,  and  then  to  linearly  interpolate  in  two  dimensions. 


22 


A  bilinear  function  of  the  form 

h(x,y)  =  a1  xy  +  a2  x  +  a3  y  +  a4 

is  constructed  so  that 

h(0,0)  =  hi;j 
h(0,l)  =  hj>j+1 

"(l-O)  1  hw,i 
hd.D  *  Vl.j+1  . 


02) 


(13) 


It  is  easily  shown  that  the  four  constants  a-|,...,a^  are  given  by 
al  =  (hij  '  hi ,j+l  ‘  hi+l ,j  +  hi+l,j+l) 
a2  =  (hi+l,j  '  hij} 


and  that 


a4  =  hij  ’ 

h(x,y)  =  h^-vxy  -  x  -  y  +  1)  +  h-+1^-(-xy  +  x)  + 

hjjtl-f-xy  *  y)  *  h1t1>j+1-(«y)  • 


(14) 


(15) 


It  should  also  be  noted  that  the  same  result  is  gotten  by  first 
linearly  interpolating  across  the  upper  and  lower  rows,  and  then 
vertically.  The  order  (rows  first  or  columns  first)  of  the  linear 
interpolations  does  not  matter. 

The  bilinear  approximation  has  the  desirable  property  that  it  is 
easily  calculated.  The  derivatives  of  h(x,y)  with  respect  to  x  and  y 
are  also  easily  found,  except  at  the  grid  points.  There  the  derivative 
is  undefined. 

Although  the  ease  of  calculation  makes  the  bilinear  approximation 
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attractive,  another  consideration  is  that  the  mathematical  analyses  that 
construct  algorithms  for  constrianed  optimization  invariably  require  in 
the  proof  of  convergence  that  the  constraint  functions  be  continuously 
differentiable.  The  bilinear  interpolation  does  not  possess  this  property, 
and  therefore  convergence  difficulties  may  arise. 

As  alternative  procedures  for  interpolation  two  other  approaches 
were  also  implemented  and  evaluated.  These  were  a  bicubic  Hermite  inter¬ 
polation  scheme  and  a  spline  based  approach. 

The  bicubic  Hermite  interpolating  polynomial  was  chosen  as 

h(x,y)  =  b]  x3y3  +  b2x2y3  +  b3xy3  +  b4y3  +  b5x3y2  +  bgx2y2  + 
b?xy2  +  b8y2  +  bgx3y  +  t>1Qx2y  +  b^xy  +  b]?y  + 

b13x3  +  b14x2  +  b15x  +  b16  (16^ 
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where  Q(x)  =  2x3  -  3x2  +  1;  Q(y)  =  2y3  -  3y2  +1;  Q(l-x)  =  2(l-x)3  -  3(l-x)2  +1 
Q(l-y)  =  2 ( 1 -y ) 3  -  3(1 -y ) 2  +  1  . 

It  is  seen  that  the  approximation  is  a  weighted  sum  of  the  four  nearest 

grid  points.  The  fact  that  the  derivatives  are  continuous  can  be  observed 

by  noting  that  h{x,y)  is  separable  and  that  dQ~^0  and  dQ_>0  as  x-*0  or  x-»l 

dx  dy 

and  as  y-^0  or  y-y  I . 

A  disadvantage  of  the  bicubic  Hermite  polynomial  interpolation  scheme 
is  that  the  derivatives  are  forced  to  be  zero  at  the  grid  points.  Although 
this  does  insure  continuity  of  the  derivative  functions,  it  is  somewhat 
contrary  to  intuition  in  that  the  terrain  at  node  points  is  made  flat. 

This  difficulty  may  be  overcome  by  using  a  larger  array  of  points  with  which 
to  calculate  an  interpolating  function. 

We  allow  the  interpolating  function  to  depend  on  the  points 
hk  m  where  k  =  i-1 ,i ,i+l ,i+2,  and  m  =  j-1 ,j ,j+l ,j+2  . 

This  is  illustrated  below  (Figure  11). 
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estimated  is 
between  x^  and  xi+1 

'  and  between  y-  and 

*j+l. 


Figure  11. 


The  figure  shows  that  the  interpolation  will  use  two  points  on  either 
side  of  the  unknown  value,  in  the  horizontal  and  the  vertical  directions 
for  a  total  of  sixteen  grid  points. 
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The  conditions  to  be  enforced  are  that  the  approximating  function 
interpolate  to  the  data  at  the  points  (xi  .y^),  (x.  »Yj+-|).  Ui+1,y.j) 

(x.j+1  ’-yj+i ) >  and  the  slope  in  the  x  and  y  directions  at  these  four  points 
satisfy  the  general  relations 


Jh  (k,m)  =  h(k+1 ,m)  -  h(k-l ,m) 
jx  2 

d_h  (k,m)  =  h(ktm+l )  -  h(k,m-l ) 

Jy  2 


k  =  i ,i+l 
m  =  J»j+1 


The  conditions  in  (20)  assure  that  the  derivative  functions  are  continuous 
and  that  the  slope  at  the  grid  points  makes  intuitive  sense. 

The  one-dimensional  version  of  this  interpolation  problem  was  termed 
convolutional  cubic  interpolation  in  a  recent  report  (Ref.  38).  The 
approximation  function  for  an  unknown  data  point  located  between  x^  and  x.+-j 

was  found  to  be  given  by 

f ( x )  =  1/2  |fi_]-(-x3  +  x2  -  x)  +  f.-(3x3  -  5x2  +1) 

+  fi+1*(-3x3  +  4x2  +x)  +  fi+2’(x3  -  *2)J. 

for  points  x>x  between  x^  and  xi+1. 

For  the  two  dimensional  case  we  define  the  functions 


a(x)  =  (-x3  +  x2  -x) 
b(x)  =  (3x3  -  5x2  +  1) 
c  ( x )  =  (-3x3  +  4x2  +  x) 
d(x)  =  (x3  -  x2) 


a(y)  =  (-y3  +  y2  -  y) 

b(y)  =  (3y3  -  5y2  +  1) 
c(y)  =  (-3y3  +  4y2  +  y) 
d(y)  =  (y3  -  y2). 


The  two  dimensional  cubic  convolution  interpolating  polynomial  is  a  weighted 
sum  of  16  data  points  and  is  most  easily  written  as  the  quadratic  form 


-1,-1 

h-l,0 

h-l,l 
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|d  ( x )] 
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It  can  be  seen  that  the  two-dimensional  convolution  cubic  requires 
approximately  four  times  as  much  work  to  evaluate  than  the  bilinear  or 
bicubic  Hermite  interpolating  polynomials.  Examples  were  run  using  each 
of  these  techniques  for  the  evaluation  of  non-grid  data  points. 

In  Section  IV  we  describe  classes  of  algorithms  which  have  been 
developed  to  solve  nonlinear  programming  problems. 
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IV,  ALGORITHMS  FOR  NONLINEAR  PROGRAMMING 


The  nonlinear  programming  problem  is  one  of  maximizing  or  minimizing 
a  functional  of  n  variables  where  the  set  of  allowable  choices  for  variables 
is  defined  by  a  set  of  equality  and  inequality  constraints.  The  equality 
constraints  state  that  the  allowable  x  values  must  lie  on  some  curve,  while 
inequality  constraints  indicate  that  the  feasible  x  values  must  lie  below 
(or  above)  some  curve.  Mathematically  the  problem  is  stated  as 

min  f(x) 

x  "  (1) 

subject  to:  h.()<)  =  0  i=l,...,me  (equality  constraints) 

g^xj^  0  j=l , . . .  ,m..  (inequality  constraints) 

It  should  be  noted  that  the  minimization  of  f(x)  is  equivalent  to  the  max¬ 
imization  of  [-f (x)]  and  that  constraints  of  the  form 
a  q(x)  3.  b, 

and  bounds  on  variables,  can  be  cast  into  the  general  form 

g^x^O  . 

The  fundamental  theoretical  result  on  the  nonlinear  programming 
problem  was  published  in  1951  by  Kuhn  and  Tucker  (ref.  39).  In  this  paper 
necessary  and  sufficient  conditions  for  constrained  extrema  were  derived. 
However,  these  conditions  did  not  directly  lead  to  iterative  solution 
methods  and  it  was  not  until  a  decade  later  that  effective  algorithms  for 
nonlinear  programming  were  developed.  In  this  section  we  describe  some  of 
the  more  successful  procedures  that  are  used  in  nonlinear  programming. 

First,  we  consider  penalty  function  methods  (refs.  13,  14,  15,  16). 

The  basic  idea  here  is  to  force  compliance  with  constraints  by  adding  a 
large  cost  term  to  an  unconstrained  functional  whenever  a  constraint  is 
violated.  The  advantage  of  this  approach  is  that  powerful  methods  for 


unconstrained  minimization  can  be  used.  For  example,  for  the  constrained 
optimization  problem  (1)  described  on  the  previous  page,  a  new  cost 
functional  could  be  formed  as: 

m.  m 

f2(x)  =  f (x)  -  rk  •  ^  {  ln£g.(x)j|+  sR-  /[hi( 

where  rR  and  sR  are  large  numbers. 

Notice  that  the  problem  is  now  unconstrained.  The  first  additional 
term  forces  x  to  stay  in  the  inequality  constrained  region,  while  the 
second  additional  term  penalizes  deviations  from  the  equality  constraints. 
The  unconstrained  value  of  f2(x.)  minimized  should  be  the  same  as  the  con¬ 
strained  value  of  f(xj  ,  provided  that  the  weights  rR  and  sR  are  chosen 
well.  It  is  often  necessary  to  solve  a  sequence  of  unconstrained  problems 
in  order  to  arrive  at  a  solution  to  the  constrained  problem.  This  is 
because  too  small  a  value  for  the  weights  will  allow  constraint  violations, 
while  too  large  a  choice  may  lead  to  numerical  difficulties  due  to  the 
insensitivity  of  f2(x)  to  changes  in  f(x).  The  Sequential  Unconstrained 
Minimization  Technique  (SUMT)  as  described  by  Fiacco  and  McCormick  in 
(Ref.  13)  has  been  among  the  most  successful  of  the  penalty  methods.  For 
an  application  such  as  the  TF/TA  problem  the  large  number  of  constraints 
would  tend  to  limit  the  usefulness  of  the  penalty  function  approach  as  it 
is  likely  that  the  sum  of  the  penalties  would  overwhelm  the  cost  that  was 
being  minimized.  On  the  other  hand,  recent  developments  in  the  parallel¬ 
ization  of  unconstrained  techniques  (ref.  17,  18,  19,  20)  may  make  the 
penalty  function  approach  a  good  one  for  problems  with  relatively  few 
constraints. 

A  highly  successful  method  (see  21  for  computational  results)  is  the 
generalized  reduced  gradient  method  developed  by  Adabie  and  Carpentier 
(ref.  22,  23,  and  applications  in  24,  25).  The  method  is  an  extension  of 


the  procedure  developed  by  Wolfe  (ref.  26)  to  the  case  of  nonlinear 
constraints.  This  method,  and  the  method  of  feasible  directions  to  be 
discussed,  provided  a  baseline  against  which  computational  improvements 
were  measured.  It  has  been  extensively  tested  and  successfully  adapted 
to  solve  problems  with  a  large  number  of  variables  and  constraints. 

The  philosophy  for  the  GRG  method  is  similar  to  that  of  linear 
programming.  Through  the  introduction  of  basis  (at  a  constraint  boundary) 
and  non-basis  variables  (within  the  boundaries)  a  simplified  reduced 
problem  is  formed.  This  reduced  problem  is  solved  and  a  direction  for 
search  is  obtained.  When  this  search  direction  is  projected  into  the 
feasible  region  (the  region  of  all  points  which  satisfy  all  constraints), 
a  direction  for  the  non-reduced  (i.  e.  the  original)  problem  is  generated. 
Several  implementations  for  the  GRG  method  have  been  developed. 

A  third  class  of  algorithms  is  the  set  of  feasible  directions  algorithms. 
This  procedure  was  first  developed  by  Zoutendijk  (ref.  27),  and  has  been 
expanded  by  Polak  (ref.  28,  29,  30,  31)  and  others  (ref.  32,  33,  34).  The 
idea  is  to  pick  a  starting  point  that  satisfies  all  the  constraints  (a 
feasible  point)  and  then  to  find  a  direction  along  which  a  small  move 
violates  no  constraints  (feasible  direction),  and  at  the  same  time  decreases 
the  cost  functional.  The  procedure  progresses  as  follows: 

min  f(x)  (3) 

subject  to  G^(x)  —0.  (we  note  that  equality  constraints 

can  be  put  in  this  form) 

The  feasible  domain  is  the  set 

F  =  ^  x  |  Gk(x)  ^  0  for  all  k }  .  (4) 

Define  a  function  PSI 

PSI(x)  =  max  jo  ;  Gk(x)j 

and  note  that  PSI(x)>0  if  x^F  and  PSI(x)  =  0  if  x  6  F. 


(5) 
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Also  define  the  "^-active  constraint  set"  as  the  set  of  all  constraints 
within  £  of  PSI(x).  That  is 

J*(x)  =  /  k  j  Gk(x)  »  PSI(x)  -6\.  (6 

The  algorithm  chooses  a  direction  d  which  is  the  weighted  sum  of  the 
(negative)  gradients  of  f(xj  and  the  ^-active  constraints  (reducing  the 
cost  vs.  moving  away  form  the  boundary).  Then  either  a  line  search  or 
sequence  of  fractional  steps  (1,  1/2,  1/4,  -  the  Armijo  rule)  is  instituted 

to  minimize  the  cost  functional  in  this  direction,  while  not  violating  any 
constraints.  Thus  a  new  feasible  point  is  found  and  the  process  is  iterated. 
Work  has  been  done  applying  the  feasible  directions  approach  to  the  TF/TA 
problem  (ref.  5),  and  while  the  procedure  does  generate  solutions  it 
appears  that  the  computation  time  may  be  too  long  for  implementation.  An 
advantage  of  the  feasible  directions  algorithm  is  that  if  the  procedure  is 
halted  before  an  optimal  point  is  found,  the  resulting  suboptimal  traj¬ 
ectory  is  still  allowable.  In  the  context  of  the  TF/TA  problem  this  would 
mean  that  the  truncated  optimization  still  leads  to  a  safe  flight  path 
that  would  be  within  the  capability  of  the  airframe  to  fly.  Some 
methods,  particularly  the  augmented  Lagrangian  methods,  do  not  have  the 
property  that  current  iterates  are  necessarily  feasible. 

The  augmented  Lagrangian  methods  are  also  called  multiplier  methods. 
Since  the  Kuhn-Tucker  conditions  are  necessary,  these  methods  attempt  to 
find  the  set  of  Lagrange  multipliers  associated  with  the  optimal  feasible 
point.  They  have  similar  structure  to  the  penalty  function  approach  but 
avoid  the  necessity  that  the  rk  and  sk  of  those  methods  approach  infinity 
in  order  to  achieve  convergence.  Fletcher  (ref.  40)  proposed  a  functional 
of  the  form 

me 

(p  (x,v,r)  =  f(x)  +  (1/2)  iL  rj[hj(x)  -  v^]2  + 
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m  +m- 


+  {1/2)  Jr  1  rjlmin[°  •  (g j(x)-vj)2]^ 


J=V1 

This  was  termed  an  exact  penalty  function  or  augmented  Lagrangian.  The 
"penalty  parameters"  r^,...,  rm  are  chosen  as  large  numbers.  The 


(7) 


'r 


me+mi 

,vm  +m  are  re^ate<1  t0  *^e  Lagrange  multipliers  and  it  is  hoped 
e  i 


that  each  v.r.  product  approximates  u*.  the  optimal  Lagrange  multiplier 

J  J  J  * 

for  constraint  j.  The  iterative  procedure  is  to  guess  initial  values  of 
the  v j  '  s ,  and  to  solve  the  unconstrained  problem 


in  (p  (x,v ,jr )  . 


mi 

x 


This  produces  a  new  x.  The  v.'s  are  then  changed  according  to  a  rule 

J 


(8) 


vj  =  Vj  -  min[  g^ ( x )  ,  y.  ] 


for  inequality  constraints 


and 


(9) 


vj  =  vj  -  ho(i> 


for  equality  constraints. 


The  unconstrained  problem  is  then  reformulated  and  solved  (equation  8 
with  the  new  .v  value).  If  the  r.'s  are  chosen  sufficiently  large  then 

J 

★ 

it  may  be  shown  that  the  x  and  v^  vectors  approach  the  optimal  x  and 

★ 

Lagrange  multipliers  u.  =  v.r.  . 

J  J  J 

Some  variations  of  the  multiplier  method  as  described  above  were 
given  by  Rockafellar  (ref.  41)  and  Pierre  and  Lowe  (ref.  42).  In  these 
formulations  the  (j)  function  (augmented  Lagrangian)  was  chosen  somewhat 
differently,  but  the  general  procedure  was  similar. 

The  final  class  of  algorithms  that  we  evaluated  was  the  quadratic 
approximation  class.  This  is  perhaps  the  most  recent  approach  and  seems 
to  be  successful.  In  1963  Wilson  (ref.  43)  proposed  that  a  constrained 
minimum  could  be  found  by  solving  a  properly  chosen  quadratic  programming 
problem  at  each  iteration.  The  cost  functional  was  a  quadratic  approximation 
to  the  Lagrangian  and  the  constraints  were  linearized  about  the  current  x. 


<  *■  ►  -*v  . . 
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The  advantage  of  such  an  approach  would  be  that  there  are  fairly 
successful  methods  for  solving  the  quadratic  programming  problem.  Wilson 
proposed  the  following  algorithm: 

Choose  a  direction  d  so  that  d  solves 

min  (1/2)  dt-[7^L(x,u)  ]-d  +  7f(x)t-d 
d  “  x  --  x  -  -  (10) 

subject  to  h.(x)  +  cryh-(x)t  d  =  0  (equality  const.) 

J  X  J 

9j(>0  +  ^xgj-(x)t  d  0  (ineq.  constr.) 

Then  =  +  <*  d 
—new  —old  — 

The  Lagrange  multipliers  for  the  quadratic  subproblem  are  used  to 

p 

define  the  new  PxL(x,u)  and  the  method  is  repeated. 

Han  (ref.  35,  page  65)  developed  a  quasi-Newton  update  for  the 
Hessian  of  the  Lagrangian  and  this  procedure  was  modified  and  implemented 
as  a  computer  program  by  Powell  (ref.  44,  45).  Another  approach  is  in  (ref.  46). 

One  of  the  goals  of  the  research  project  was  to  identify  computational 
procedures  that  would  be  applicable  to  large  constrained  optimization 
problems.  For  example,  line  searches,  the  solution  of  reduced  sets  of 
nonlinear  equations,  quadratic  programs,  etc.  are  inherent  steps  in  the 
overall  optimization  process  described  by  each  of  the  classes  of  algorithms. 
Particular  methods  to  accomplish  such  steps  change  the  performance  of  the 
general  algorithm.  Some  work  in  parallelization  has  been  reported  by 
Mukai  (in  ref.  17)  and  by  Chazen  and  Miranker  (ref.  18),  while  the 
quasi -Newton  concept  developed  by  Broyden  (ref.  36)  has  been  applied  to 
constrained  optimization  by  Han  (ref.  35,  page  65). 

Computer  programs  which  implement  particular  versions  of  the  general¬ 
ized  reduced  gradient,  augmented  Lagrangian,  and  quadratic  approximation 
methods  were  obtained  and  slightly  modified.  These  programs  are  described 
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in  Section  V  of  this  report.  An  attempt  to  develop  an  implementation 
of  the  feasible  direction  method  described  by  Polak  (ref.  31,  33)  proved 
unsuccessful.  No  attempt  was  made  to  test  penalty  function  methods  since 
these  are  now  considered  noncompetitive  (ref.  21). 
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V.  COMPUTER  PROGRAMS  FOR  CONSTRAINED  OPTIMIZATION 

Several  programs  for  solving  the  general  nonlinear  programming 
problem  were  obtained.  One  generalized  reduced  gradient  procedure,  two 
multiplier  methods  (  augmented  Lagrangian),  and  one  quadratic  approxi¬ 
mation  method  were  evaluated  in  a  series  of  tests.  In  this  section  we 
describe  the  attributes  and  the  usage  of  these  programs. 

The  first  procedure  we  implemented  was  a  program  that  performed  the 
generalized  reduced  gradient  algorithm.  The  program  is  named  GRG  and  was 
developed  by  Lasdon  (ref.  24),  and  has  undergone  several  revisions  since 
that  time.  The  program  is  structured  as  a  stand  alone  subroutine,  and 
all  problem  data  (constraint  bounds,  variable  bounds,  etc.  )  is  read  in 
from  the  card  reader  stream.  A  user  supplied  subroutine  CALCFG  is 
required  and  contains  the  calculations  of  the  objective  function  and  the 
constraints.  The  user  may  specify  whether  derivatives  of  objective  and 
constraints  are  to  be  supplied  analytically  through  a  user  supplied 
subroutine  PARSH,  or  whether  these  are  to  be  calculated  by  finite  difference 
approximation.  The  program  consists  of  seventeen  subroutines  and  the 
FORTRAN  code  is  approximately  262b  statements  in  length. 

The  actual  form  of  the  problem  solved  by  GRG  is: 
min  gM,..(x) 

x  M+1  (1) 

subject  to  g.(x)  =  0  i  =  l,...,  NEQ  (equal,  constr. 

0  ^g^x)  £  UB(N+i)  i=NEQ+l,...,M 

(ineq.  Constr.) 

LB(i)i  x.  ^  UB(i )  i=l , . . .  ,N 

That  is:  There  are  N  variables,  NEQ  equality  constraints,  M-NEQ  inequality 
constraints,  and  each  variable  may  have  lower  and  upper  bounds  associated 


with  it. 
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The  program  adds  slack  variables  x^+1  , . ,XN+M  to  the  set 

of  "natural  variables"  (x-j .  ,x^) ,  and  produces  the  problem: 

7  ^  ,2, 

subject  to  g^  (x)  -  x^+i  =0  i=l,...,M 

LB(i)  6  x.  *=  UB(i)  i  =1 , . . .  ,N+M 

where 

LB ( i )  =  UB(i)  =  0  i=Nfl , . . . ,N+NEQ 
LB( i )  =  0  i=N+NEQ+l . . ,  N+M. 

The  procedure  works  as  follows.  Suppose  the  current  value  of  x 
satisfies  all  constraints,  and  that  a  number  NB  of  these  are  binding 
constraints.  The  program  uses  the  NB  binding  constraint  equations  to 
solve  for  NB  of  the  natural  varaibles  in  terms  of  the  other  N-NB  natural 
variables  and  the  NB  slack  variables  that  are  associated  with  the  binding 
constraints.  These  N  variables  (N-NB  +  NB)  are  called  nonbasic  variables. 
Now,  let  w  be  the  vector  of  basic  variables  (length  NB)  and  y  be  the 
vector  of  nonbasic  variables  (length  N).  The  vector  of  binding  constraint 
functions  can  then  be  denoted 

£(w,y  )  =  0  .  (£  is  of  length  NB)  (3) 

The  set  of  basic  variables  can  be  chosen  so  that  the  matrix  B  (NB  x  NB) 

B  4  (^./dWj)  (4) 

is  nonsingular  at  the  current  value  of  x. 

Using  (4),  the  value  w  in  (3)  may  be  solved  for  in  terms  of  y.  The 
objective  function  is  therefore  a  function  of  y  only  (N  variables) 
and  this  reduces  the  original  problem  (N+M  variables)  to  a  simpler  problem 
min  gM+1(w  {y),i)  =  F(y) 

y  (5) 

subject  to  l_b  £  y  &  ub  . 

The  original  problem  (1)  is  solved  by  generating  and  solving  a  sequence 


36 


of  reduced  problems.  Notice  that  the  reduced  problems  have  no  nonlinear 
constraints  and  so  can  be  solved  by  a  gradient  method.  The  search 
directions  for  the  reduced  problem  are  determined  by  the  BFGS  update  and 
a  quadratic  polynomial  fit  is  used  for  one-dimensional  line  searches. 

Once  the  search  direction  is  found,  a  line  search  is  performed  to 
minimize  F(y^  +^d  )  and  this  requires  the  solution  of  the 

nonlinear  equations 

£  (w  ,  i  +  *  d  )  =  0.  (6) 

It  is  possible  that  some  of  the  basis  variables  will  now  violate  their 
bounds.  If  this  is  so,  a  new  set  of  basis  variables  will  be  picked  and 
a  new  reduced  problem  created. 

For  the  6RG  program,  the  initial  guess  for  the  solution  vector  x 
need  not  be  feasible.  The  program  will  create  a  subproblem  that  finds 
a  feasible  point  if  one  is  not  provided.  At  the  conclusion  of  each 
iteration  of  the  main  procedure  the  method  guarantees  that  the  current 
value  of  x  is  feasible,  and  that  the  objective  function  has  been  decreased. 

The  second  program  that  was  implemented  was  LPNLP,  as  described 
in  reference  (42).  This  procedure  is  an  augmented  Lagrangian  method 
and  the  FORTRAN  code  consists  of  approximately  1200  statements.  The 
program  is  organized  as  a  subroutine  to  be  called  by  a  main  program.  The 
main  program  dimensions  all  arrays  and  workspaces.  Two  user  supplied 
subroutines  must  be  provided.  Subroutine  FXNS  provides  the  objective 
function  and  the  constraint  equations.  Subroutine  GRAD  provides  the 
derivatives  of  the  objective  and  the  constraints.  No  mechanism  for 
the  numerical  calculation  of  derivatives  is  provided,  and  due  to  the 
intractability  of  these  calculations  for  the  test  problem,  a  finite 
differencing  scheme  was  used  to  produce  the  data  needed  in  subroutine 
GRAD.  All  constraint  bounds  and  bounds  on  variables  were  read  in  on 
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the  card  reader  stream.  Upper  and  lower  bounds  on  all  variables  were 
explicitly  treated.  The  actual  form  of  the  problem  solved  by  LPNLP  was: 

max  f(x) 
x 


subject  to 

hi(x)  =  ai 

i=l . .  NE 

9j(x)  *  bj 

j=l,....,  NI 

and 

ck  -xk  -dk 

k=l, _ ,n  . 

The  procedure  used  is  to  form  an  augmented  Lagrangian  of  the  form 

NE  NI 

L  (x,  £,w)  =  f(x)  +  ^  xi  (ai “h-  (x))  +  £  ^-(b.-g.(x)) 

i=l  j=l  J  J  J 

NE  os' 

-  w,  ^  [a.-h. (x)r  -  w?  [b.-g-(x)]' 

1  i=l  1  1  L  j  £  C  J  J  _ 

3 

-  w3  ^  [b.-g.(x)]2  (8) 

j£CK  J  J 


c -  LL».-y.| 

0  i/r  J  J 


where  /  /  .  .• 

ca  -  /j  |A  >  0} 

Cb  =  /j  |  4j  =0  and  g-j(x)^bj}. 

The  values  w-j ,  w^,  and  w3  are  penalty  weights,  while  the  <^'s  and  ^'s 
are  Lagrange  multipliers.  Thus  the  first  three  terms  of  the  expression 
correspond  to  the  Lagrangian  for  the  constrained  problem  and  the  next 
three  terms  are  penalties  on  constraint  violations. 

The  method  picks  initial  «s  ^9°  and  we,  and  solves  the  unconstrained 
problem 

max  L  (x  ,  pc*,  w*).  ,q, 

x  v 

Let  be  the  solution  to  the  problem  (9).  This  vector  will  normally  not 

satisfy  all  the  constraints.  At  x^  the  gradient  with  respect  to  x 

^x  La^’  »  (1C 

whereas  the  Kuhn-Tucker  conditions  for  the  original  constrained  problem 


require  that 


V»L(x,  *,<£)=  0.  (H) 

The  procedure  reassigns  *.  and  B  so  that 

L(x]_  ,  S')  =  ^.L^x1,  ±  ,4  *  w  )  *  0  .  (12) 

At  this  point,  the  weights  w  may  be  increased.  The  new  Lfl(x,  *.  ,  B  ,  w  ) 
is  then  maximized  with  respect  to  x  and  the  update  process  for  and  B 
(the  Lagrange  multipliers)  is  repeated. 

The  program  provides  several  options  that  may  be  employed  in  the 
generation  of  the  search  direction  for  the  unconstrained  maximization. 

The  general  procedure  for  this  is  the  Davidon-Fletcher-Powell  update; 
however,  one  may  specify  a  self-scaling  variable  metric  method  (ref.  42 
page  405).  Additional  flexibility  is  provided  in  that  a  mode  is  possible 
in  which  the  DTP  direction  is  reset  to  the  gradient  direction  every  N 
searches . 

One  disadvantage  of  this  and  the  other  augmented  Lagrangian  methods 
is  the  problem  of  constraint  breakthrough.  The  current  value  of  x  at 
the  end  of  an  iteration  is  not  explicitly  forced  to  be  feasible.  However, 
by  allowing  the  penalty  weights  to  be  increased  as  the  iterations  progress 
one  may  limit  this  problem. 


A  third  program  we  implemented  was  VF01AD.  The  program  is  a 
member  of  the  Harwell  Subroutine  Library  and  was  written  by  R.  Fletcher 
as  described  in  (ref.  40).  It  consists  of  approximately  940  statements 
of  FORTRAN  code.  The  program  is  an  augmented  Lagrangian  method,  and  the 
problem  solved  is 


min  f(x) 
x 


(13) 


subject  to  h^ (x)  =  0  i=l,...,NEQ 

g.(x)£0  j=NEQ+l . .  M. 

J 
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In  this  procedure  the  augmented  Lagrangian  is  given  by 

„  NEQ  2 

9  )  =  f(x)  +  J  ^  (h^  (x)  -  9^  y  + 

M  ? 

+  i  ZL  min[  0  ,  (gi (x)  -  9.)  ]  . 

j=NEQ+l  J  J  -  J 

Each  iteration  involves  the  unconstrained  minimization  of  this  functional 
for  a  fixed  9  and  ^  .  At  the  end  of  each  iteration  these  are  changed  so 
that  the  sequence  of  unconstrained  minima  tend  to  the  solution  of  the 
constrained  problem  (13).  At  solution,  the  products  0^.  are  the 
Lagrange  multipliers  corresponding  to  each  of  the  constraints  i=l,...,M. 
Convergence  is  guaranteed  (with  exact  arithmetic)  and  the  method  can  be 
expected  to  converge  at  a  second  order  rate  (ref.  40). 

An  initial  guess  for  the  solution  vector  must  be  provided,  but  this 
need  not  be  a  feasible  point.  The  user  must  also  supply  a  subroutine 
which  calculates  the  objective  functional,  the  constraints,  and  derivatives 
of  the  objective  and  constraints.  Linear  constraints  are  handled  in 
an  efficient  manner,  and  prior  information  about  the  form  of  the  Hessian 
matrix  or  the  Lagrange  multipliers  can  be  utilized.  The  unconstrained 
minimizations  are  performed  by  a  subroutine  VAQ9AD  which  is  an  efficient 
quasi -Newton  method  that  avoids  line  searches. 

Finally,  we  obtained  the  quadratic  approximation  method  of  Wilson, 

Han  and  Powell  (refs.  35  page  65;  43,  44).  This  program  (VF02AD)  was  written 
by  Powell  and  is  part  of  the  Harwell  Subroutine  Library.  It  consists 
of  approximately  1183  lines  of  FORTRAN.  The  problem  solved  is  of  form 

min  f(x) 
x 

h .  (x)  =  0  i=l . NEQ 

g.(x)  £  0  i=l . M-NEQ  . 


subject  to 


(15) 


The  method  of  solution  is  iterative  where  each  iteration  minimizes  a 
quadratic  approximation  to  the  Lagrangian  subject  to  linear  approximations 
to  the  (nonlinear)  constraints.  The  constraints  in  the  quadratic  programm¬ 
ing  subproblem  are  modified  to  avoid  infeasibility.  At  each  iteration  the 
subproblem  is  dependent  on  the  current  x  and  is 
min  d^^  U*  ,  u  )*d  +  C^flx)1  d 

subject  to  h.(x)  +  P”  h -(x)^-  d  =  0  (equality  constr. ) 

J  X  J 

g.(x)  +  (7  g.  (x^-  d  ■i  0  (ineq.  constr.) 

1  X  J 

The  direction  d  is  found  by  solving  the  quadratic  program  (10)  and  the 

new  value  of  is  chosen  as 

x  =  x  ,  .  +  «.  d  . 

-new  —old  — 

The  stepsize  is  determined  by  a  quadratic  interpolation  of  a  penalty 

2 

function,  and  the  program  automatically  estimates  17 ^  L(x,u)  by  an 
extension  of  the  variable  metric  method  of  unconstrained  minimization. 

The  quadratic  program  that  forms  the  subproblem  is  solved  by  the  routine 
VE02AD. 

The  quadratic  approximation  method  as  exhibited  here  has  an  advantage 
over  augmented  Lagrangian  methods  in  that  constraints  are  explicitly 
considered  in  the  calculation  of  the  search  direction.  Thus,  feasibility 
of  the  current  iterate  is  somewhat  assured. 

It  has  been  stated  (ref.  44)  that  although  more  work  is  required  to 
solve  the  quadratic  subproblems  of  this  method  than  the  unconstrained  min¬ 
imizations  of  the  augmented  Lagrangian  methods,  the  total  number  of 
function  and  gradient  calculations  should  be  much  lower.  This  may  be  an 
important  consideration  when  the  objective  and  constraint  calculation  is 
costly. 

The  user  of  VF02AD  is  required  to  provide  an  initial  guess  for  the 


solution  vector.  This  guess  need  not  be  a  feasible  point.  He  must 
also  supply  code  to  evaluate  the  objective  functional,  constraint 
functions,  and  the  derivatives  of  these. 

Each  of  the  above  programs  was  compiled  on  a  Digital  Equipment 
Corporation  VAX  11/780  computer.  Three  of  the  programs  (LPNLP,  VF01AD, 
and  VF02AD)  were  modified  so  that  derivatives  could  be  calculated  by 
finite  difference  approximations.  The  program  GRG  already  had  this 
facility  built  into  it.  Several  of  the  programs  had  machine  dependent 
parameters  or  assembly  laneuage  portions  which  had  to  be  altered  before 
they  could  be  run  on  the  VAX  system. 

We  attempted  to  write  a  program  to  implement  the  feasible  directions 
algorithm  of  Polak  (refs.  30,  31)  but  were  unable  to  make  it  reliable.  This 
was  undoubtedly  the  fault  of  this  programmer,  not  of  the  algorithm.  Also, 
for  completeness  we  mention  two  other  programs  we  had  hoped  to  test.  The 
first  of  these  is  a  quadratic  approximation  technique  that  is  different 
from  that  of  VF02AD.  The  program  is  called  0PRQP/XR0P  and  was  written 
by  M.  C.  Bartholomew-Biags  (ref.  46).  Unfortunately  we  were  unable  to 
obtain  a  copy  of  this  from  the  author  in  time.  Second  was  a  program 
called  MINOS  which  was  furnished  by  the  Stanford  University  Center  for 
Information  Technology.  This  9600  line  program  is  designed  to  handle 
large,  sparse  nonlinear  programming  problems.  The  algorithm  is  a  projected 
augmented  Lagrangian  method  which  solves  a  sequence  of  linearly  constrained 
subproblems  by  a  reduced  gradient  technique.  We  received  this  package 
too  late  to  adapt  our  test  problem  to  the  rather  intricate  input  format 
required.  The  presence  of  several  undocumented  machine  dependent  parameters 
also  hindered  the  use  of  this  program  in  our  study. 

In  Section  VI  we  describe  our  test  problem  and  the  computational 
results  for  each  of  the  subroutines  we  studied. 
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VI.  COMPUTATIONAL  RESULTS 

In  this  section  we  describe  a  series  of  tests  that  were  used  to: 

1 )  Check  the  correctness  of  the  computer  code  for  the  four  programs 
that  were  obtained 

2)  Demonstrate  that  the  inclusion  of  a  finite  differencing  scheme 
for  gradient  calculations  was  correct 

3)  Assess  the  relative  performance  of  the  four  programs  on  small  and 
medium  sized  nonlinear  programming  problems  which  had  known  solutions 

4)  Evaluate  the  feasibility  of  attempting  to  solve  the  path  generation 
problem  of  terrain  fol lowing/terrain  avoidance  flight  path  generation 
in  real  time 

and 

5)  To  modify  the  formulation  of  the  path  generation  problem  in  an 
attempt  to  improve  the  speed  of  solution. 

As  stated  in  Section  V,  four  computer  programs  were  used  in  the  stuoy. 
In  order  to  make  sure  that  these  were  correctly  implemented  on  our  system, 
and  to  become  familiar  with  their  use  (choice  of  parameters,  input 
formats,  etc.)  we  first  attempted  to  find  constrained  minima  of  several 
standard  problems  whose  solutions  are  known.  A  typical  example  was 

p  p  p 

min  f ( x)  =  x^  +  x2  +  x^x2  -1 4x1  -  16x2  +  (x3  -  ID)  + 

4(x4  -  5)2  +  ( x 5  -  3)2  +  2(x6  -  l)2  +  5x72  + 

7(x8  -  ll)2  +  2(xg  -  10)2  +  (x1Q  -  7)2  +  45 

subject  to  the  set  of  constraints 

105  -  4x-j  -  5x2  +  3x7  -  9xg  1?  0 
-10x^  +  8x2  +  17x7  -  2xg  ^  0 


8x-j  - 

2*2  - 

5xg  + 

2xiq  +  12  0 

-3(x1 

-  2)2 

-  4(x2 

-  3)2  -  2x32  +  7x4  +  120 

>  0 

-5xi2 

-  8x2 

(x3 

-  6)2  +  2x4  +  40  0 

-0.5(x]  -  8 

)2  -  2( 

x2  -  4)2  -  3x^2  +  Xg  +  30 

<-  0 

2 

'X1  ' 

■  2(x2 

-2)2 

+  2x-jX2  -  14x5  +  6xg  0 

3x-j  - 

6x£  - 

12(x9 

-  8)2  +  7x1q  V  0  . 

This  problem  has  ten  variables  and  eight  inequality  constraints.  The 
objective  function  to  be  minimized,  as  well  as  four  of  the  constraints, 
are  nonlinear.  A  feasible  point  given  by 

(*1* . . x i o )  =  (2,3,5,5,1,2,7,3,6,10) 

was  chosen  as  an  initial  guess  for  the  solution  vector.  At  this  value 
of  x  the  objective  function  has  value  f(xj  =  753. 

The  exact  solution  to  the  problem  has  a  minimum  objective  function 
f(x* )  =  24.3062091  . 

All  of  our  programs  converged  to  the  solution  successfully.  TaMe  1 
compares  solution  times,  function  evaluations,  etc. 

TABLE  1. 

Relative  Success  of  the  Four  Programs 
in  Solving  (1) 


Program 

name 

Outer  Iterations 
performed 

Total  Function 
evaluations 

Execution  Time 
(double  precision) 

GRG 

18 

367(220) 

4.35  sec 

LPNLP 

43 

583(430) 

6.6  sec 

VF01AD 

5 

792(720) 

4.32  sec 

VF02AD 

11 

154(140) 

10.09  sec 

Under  Total  Function  Evaluations  the  two  numbers  refer  to  the  actual 
number  of  times  the  objective  function  and  constraints  were  calculated,  and 
the  amount  of  that  number  that  were  part  of  gradient  calculations  by  finite 
differences.  Therefore,  the  number  of  gradient  calls  were  22,  43,  72,  and  14. 
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From  the  Table  the  most  striking  fact  is  that  VF02AD  took  substantially 
longer  than  the  others,  but  used  far  fewer  function  evaluations.  In  a 
situation  where  the  evaluation  of  the  function  is  a  time  consuming  operation 
this  program  may  have  an  advantage.  The  author  (ref.  44)  states  that  this 
program  should  be  considered  a  preliminary  version  since  the  algorithm  for 
solving  the  quadratic  subproblems  has  not  been  optimized  for  this  particular 
application. 

We  also  notice  that  the  generalized  reduced  gradient  algorithm  (GRG) 
had  significantly  fewer  function  evaluations  than  did  VF01AD  (  an  augmented 
Lagrangian  method),  and  these  were  similar  in  execution  time. 

The  relative  rankings  among  the  four  routines  as  shown  in  the  table 
were  fairly  consistent,  with  GRG  usually  being  the  fastest  and  VF02AD  the 
slowest.  The  routine  LPNLP  seemed  to  be  sensitive  to  the  particular  problem, 
and  for  some  cases  it  had  not  reached  a  solution  in  the  allotted  maximum 
number  of  iterations. 


The  next  step  in  our  investigation  was  to  design  a  test  problem  of  a 
form  that  matched  the  TF/TA  path  generation  problem  in  its  format.  The 
continuous  time  version  of  the  problem  was 


min  j  f  j  pz(t)2  +  W-[py(t)  -  yd(t)]2j  dt 

p(t),nz(t)  0 

subject  to  the  state  equations  (derived  for  a  point  aircraft) 


Px  =  V  cos  >  cos  P 
py  =  V  cos  1  sin  f’ 
Pz  =  V  sin  ^ 


three  dimensional  position 


H*  =  -(nz  sin^)/(V  cos  ^ ) 

4 

$  =  (n  cos  ft  -  g  cos  Y)  /  V 

4  ^ 

^  =  P(t) 

V  =  - g  sin  Y 


heading 

flight  path  angle 
roll  rate 
velocity 


(2) 


At  each  time  t  the  following  bounds  must  hold 


p7(t)  terrain  height  plus  a  safety  factor 


nnn  K  '  max 


max 


max 


n7  .  ^  nz(t)  “  n2 


flight  path  angle  bounds 

bank  angle  bounds 

roll  acceleration  bounds 
load  factor  constraints 


min 


max 


•  •*  *4 

nz  -  nz(t) -  pitch  jerk  constraints 

min  max 

The  control  variables  for  the  problem  are  the  load  factor  (nz(t))  and  the 
roll  rate  (p(t)).  Tne  number  W  is  a  weighting  factor  between  deviations 
from  a  desired  ground  path  y^(t),  and  altitude  minimization.  A  large 
value  of  W  forces  the  path  to  be  close  to  y^,  while  a  small  weighting 
factor  allows  deviations  if  a  lower  altitude  path  can  be  found.  The 
weighting  factor  could  also  be  a  function  of  time  or  position. 

In  order  to  put  the  problem  above  into  a  framework  that  is  solvable 
by  mathematical  programming  methods  we  discretize  the  problem  in  time. 


Let  the  interval  t=0  to  t=t^  be  broken  into  a  set  of  N  steps  of 
length  h.  That  is:tf  =  N-h.  Now,  if  the  control  variables  n z ( t )  and 
p(t)  are  considered  to  be  constants  over  intervals  of  the  form 
(  k  h  ,  k+1  h  )  i.  e.  one  time  step, 
the  Euler  integration  step  may  be  used  to  turn  the  original  differential 
equations  into  a  set  of  difference  equations. 

Euler's  method  is 

if  x  =  f(x,t)  then  x([k+l]  h)  ~x(k-h)  +  h*f(x(k*h) ,  k- h) 
provided  h  is  small.  Starting  at  k=0,  an  approximation  to  the  solution  of 
the  differential  equation  is  given  by  the  sequence  x(j-h),j=0 . ,N  . 


There  are  sets  of  bounds  on  f  and  another  set  on  nz  .  We 
state  these  in  terms  of  the  variables  by  using  finite  differences. 

n2(k)  ~  (nz(k+2)  -  2  nz(k+l)  +  n2(k)  )/  h2  k=l,...  ,N-2 

(5) 

and  p(k)  "  (p(k+l)  -  p(k)  )/  h  k=l - - ,N-1 

Finally,  there  is  the  altitude  constraint.  This  is  given  by 

p  (k)  -  [terrain(  p  (k),p  (k)  )  +  safety  factor]  £  0  k=l,...,N  . 

2  x  y 

The  cost  fucntional  is  approximated  by  the  sum 

N 

^  |pz(i)2  +  W-[Py(i)  -yd(i)  fj-  (6) 

It  may  be  seen  that  as  stated  the  problem  is  one  of  minimizing 

(12  N  +  2)  variables,  with  7  N  equality  constraints,  5  N  -  6  inequality 
constraints,  and  6  N  variables  with  bounds. 

The  values  of  the  bounds  were  chosen  as  typical  values  and  are 
shown  below: 

tf  =  flight  path  angle  =  between  -.2rad  and  ,5rad 
Pz  =  altitude  =  greater  than  60  feet  above  sea  level 
py  =  constrained  to  be  within  600  feet  of  yd 
<P  =  roll  angle  =  -lrad  to  +lrad 
nz  =  load  factor  =  between  0,3g  and  3  g 

nz  =  pitch  jerk  =  between  -1.5g  and  3g 

p  =  roll  rate  =  -lrad/sec  to  +lrad/sec 

The  initial  velocity  was  chosen  as  500ft./sec,  and  terrain  clearance^- 50  ft. 

It  is  usually  desirable  to  scale  an  optimization  problem  so  that  the 
objective  functional,  constraints  and  variables  take  on  values  of  approx¬ 
imately  the  same  order  of  magnitude.  Scaling  to  an  approximate  magnitude 
of  1  leads  to  improved  numerical  properties. 
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Let  as  choose  scale  factors 


Px ( new ) 

=  px/iooo 

Py(new) 

=  Py/1000 

Pz(new) 

*  Pz/500 

y  (new) 

=  5  r 

(new) 

=  5  y 

V  (new) 

=  V/500 

nz(new) 

=  nz/g  =  nz/32.2 

p(new) 

=  5  p 

* 

nz(new) 

=  nz/32.2 

nz(new) 

=  nz/32.2 

p  (new) 

=  5  p 

ft  (new) 

=  5  <P 

This  has  the  effect  of  changing  the  equality  constraints  and  the  bounds 
on  variables. 

For  the  terrain  function  we  first  chose  a  smooth  function  of  the 

form 

terrain(x,y)  =  800  exp(-vali  )  exp(-va!2^) 

where 

vail  =  ( x-2500 ) /1000  ;  va!2  =  (y-1000)/250  . 

This  corresponds  to  an  800  foot  hill  centered  at  (x,y)  =  (2500,1000). 

Since  straight  and  level  flight  would  consist  of  the  sequence  of 

controls  nz(t)  =  1  and  p(t)  =  0,  the  initial  guess  provided  by  evaluating 
the  state  trajectory  due  to  such  a  sequence  of  controls  and  initial 
altitude  850ft  certainly  is  a  feasible  point  in  the  space  of  variables. 

A  five  segment  problem  was  attempted.  We  found  that  even  for 

this  small  problem  none  of  the  programs  converged  to  a  solution  in  a 

reasonable  amount  of  time  (20mins.  of  computer  execution).  This  forced  us 


to  look  again  at  the  problem  we  had  stated.  We  noted  that  although 
the  state  variables  of  the  airplane  were  considered  free  variables 
and  then  were  subject  to  a  set  of  equality  constraints  that  arise  from 
the  state  equations,  an  alternative  point  of  view  would  be  that  only 
the  control  variables  (n  and  p)  are  free,  and  that  these  are  the  variables 
to  be  optimized  over.  The  state  variables  are  not  free  but  rather  are 
forced  once  the  control  sequence  and  initial  conditions  are  given.  That 
is:  the  objective  functional  and  the  equality  constraints  are  implicit 
functions  of  nz  and  p.  The  state  variables  and  state  equations  serve 
only  as  "black  boxes"  in  the  generation  of  the  inequality  constraints 
and  boundedness  conditions.  The  output  of  the  diagram  below 

State  |  Constraint  constraints 

Equations  state  I _ Functions  g(x(u)) 

variables 

is  also  seen  as  a  function  of  control  sequence  only. 

u  j  the  whole  g2(u) 

process  * 

Having  made  that  observation  we  reformulated  the  problem  as  one  of 
selecting  2  N  control  variables,  subject  to  12  N  -  6  inequality  constraints, 
no  equality  constraints,  and  such  that  there  are  N  upper  and  N  lower 
bounds  on  the  variables. 

The  scaled  and  reformulated  problem  is  shown  on  the  following 
page  (Fig.  12). 

With  the  reformulated  problem  we  obtained  convergence  in  a  five 
step  problem  with  each  of  our  four  programs.  Since  the  five  step 
problem  is  unrealistic  and  uninteresting,  most  of  our  tests  were  done 
on  problems  involving  10  or  more  steps. 
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SUBROUTINE  FXNS(X,F,FE,FI) 

IMPLICIT  REAL*8(A-H ,0-2) 

DIMENSION  X( 1 ) ,FE( 1) ,FI( 1 ) ,YDES(20> ,XX(201 ) 

K-10 

T-1.0 

GRAV=.322 

K=//0F  TIME  STEPS,  T=TIM£  STEP  SIZE 
DO  22  1=1, K  . 

ydes(I)=1.  I  the  desired  y  path  is  y-IOOOft 

vo=i. 


generate  first  step  of  state 
equations 


xo=o.o  ! 

Y0=l.l0  S  initial  values 

20=1 . 7  v=500,  y=1100,  z=850 

WEIGHT=0 .01 
XX( 1 )=(XO+T*(VO/2 . ) ) 

XX(K+1)=Y0 

XX(2*K+1 )=Z0  .  , .  , 

xx( 3*K+l )=0 .0  generate  first 

XX(4*K+l )=(T/VO)*GRAV*(X( 1 )-l . )  J  e1L 

XX(5*K+1)=T*X(K+1)  - 

XX(6*K+1)=V0 
KM1=K-1 
SUM=0 .0 
DO  1  1=1, KM1 
VEL0C=XX(6*K+I) 

PSI=XX(3*K+I)*.2 
GAMMA=XX( 4*K+I )* . 2 

phi=xx(5*k+I)*.2  calculate 

ENZ=X( 1+1 )  sta 

PEE=X(K+I+1) 

XX(  1+1 )=(XX( I)+T*( VELOC/2 . )*DCOS( GAMMA )*DCOS( PSI) ) 

XX( K+I+ 1 ) = ( XX( K+I )+T* ( VELOC/ 2 . ) *DCOS( GAMMA ) * 

&DSIN( PSI ) ) 

XX(2*K+I+1 )=(XX( 2*K+I)+T*VEL0C*DSIN( GAMMA) ) 
XX(3*K+I+1)=(XX(3*K+I)-(T*ENZ*GRAV* 

&DSIN( PHI )/ (VELOC*DCOS( GAMMA) ) ) ) 

XX(4*K+I+l)=(XX(4*K+I)+(T*GRAV/VELOC)*(ENZ* 

&DC0S ( PHI ) -DCOS ( GAMMA) ) ) 

XX( 5*K+I+1 )=( XX( 5*K+I)+T*PEE) 
XX(6*K+I+l)=(XX(6*K+I)-T*GRAV*DSIN(GAMMA)/5.) 
SUM=SUM+(XX(2*K+I)**2)+WE1GHT*4.*(XX(K+I)-YDES(1))**2 
CONTINUE 

F=SUM+XX(2*K+K)**2+WEIGHT*4.*<XX(K+K)-YDES(K))**2 
F=-F 

DO  19  1=1, K 
FI(I)=-(XX(K+I)-.6) 

FI(K+I)»-(XX(2*K+I)-(60./500.))  / 

FI(2*K+I)=-(XX(4*K+I)+1.)  I  con 

FI(3*K+I)=-(XX(5*K+I)+5.)  \ 

IF(I.EQ.K)  GO  TO  19  \ 

FI(4*K+I)=-((X(K+I+l)-X(K+I))/T+5.) 

IF(I.EQ.KMl)  GO  TO  19 

FI( 5*K+I-1 )=- ( ( X( 1+2 )-2 . *X( 1+1 )+X( I) )/(T**2)+l .5) 

CONTINUE  / 


calculate  the  sequence  of 
state  variables 


objective 

function 


constraints  are 
functions  of  the 
state  variables 


Fig.  12(a) 
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IVAL=6*K-3 

ClhHHhHhHHtlHHHHtTRlS  IS  THF.  INEQUALITY  CONSTRAINT  FOR  ALTITUDE#////////// 
DO  4  11=1, K 
XII=XX(II)*1000. 

XKP I I=XX( K+ 1 1 ) *  1 000 . 

CALL  TERRAIN(XII , XKP I I ,TERR) 

F I ( IVAL+II)=-(XX( 2*K+II )-( (TERR+50 . )/ 500 . ) ) 

4  CONTINUE 

C //#//?///// ////// //END  OF  ALTITUDE  CONSTRAINT////////////////////////// 

IVL2=6*K-3 

IVAL=7*K-3 

JJ=IVAL 

ISK1=K 

ISK2=2*K 
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DO  99  11=1, IVL2 

IF( II .GT. ISK1 .AND . II .LE . ISK2 )  GO  TO  99 
JJ=JJ+1 


FI(JJ)=-FI(II) 

CONTINUE 

RETURN 

END 


fig  12  (cont) 


Gradients  of  the  cost  function  and  the  constraint  functions  were 


calculated  by  finite  difference  approximation.  A  program  which  illustrates 


this  is  shown  in  Fig.  13.  The  derivative  of  a  function  Cj(x)  with  respect 


to  variable  xi  is  approximated  by 


20 


35 

3 


j C j (X)  ^  Cj(x  +  h-j ) 


SUBROUTINE  VF01BD(N,M,X) 

IMPLICIT  REAL*8(A-H,0-Z) 
DIMENSION  X( 1 ) ,C2( 134) 

COMMON  /VF01CD/F 
COMMON  /VF01DD/G(50) 

COMMON  /VF01ED/C(402) 

COMMON  /VF01FD/GC(20,134) 

K=N/2 

CALL  V2CONST(K,N,M,X,C,F) 

DO  3  1=1, N 
TEMP=X( I ) 

DELX* 1 .0D-05 

IF( DABS ( TEMP ).LT.1.0D0)  GO  TO  20 
DELX=DELX*DABS ( TEMP ) 

X( I )=TEMP+DELX 
CALL  V2C0NST(K,N,M,X,C2,F2) 
G(I)=(F2-F)/DELX 
DO  35  J-1,M 

GC(I,J)-(C2(J)-C(J))/DELX  > 

X( I > “TEMP 
RETURN 
END 


where  h-  is  zero 
except  for  a  small 
value  in  the  i  place 


-  evaluate  the  objective  and 

the  constraints.  V2C0NST 
is  similar  to  that  of  Fig.  12 


-  change  each  element  of  the 

variable  vector  by  a 
small  amount 

-  the  gradient  is  calculated 

-  the  variable  is  returned  to 

its  original  value 


Fig,  13  -  Finite  Differences  for 
Gradient  Calculations 

-  I-  . 


For  the  terrain  constraint 


Pz  -  terr(px  ,  p  )  -  50  ^  0 
we  evaluate  the  gradient  by  the  chain  rule, 


<j  P 


Jxi 


l  -  d  terr(px,p  ) 


dxi 


■f 


=  _/ dterr(px,py) 


d  Pv 


j  px  +  dterr(px,py)  ^ 
dXi  dPy 


+  dPZ 

The  partial  derivatives  of  the  terrain  function  are  found  by  analytically 
differentiating  the  terrain  function  (i.  e.  calculating  the  slope  at  the 
point  (PxiPy)  )•  The  other  partial  derivatives  are  found  by  the  finite 
difference  method  as  in  Figure  13. 


With  the  new  formulation  for  the  problem  we  were  able  to  obtain 
convergence  for  each  of  our  four  test  programs.  Typical  results  appear 
in  Table  2,  which  treats  a  10  step  problem.  For  this,  the  number  of 
variables  is  20,  the  number  of  inequality  constraints  is  114,  and  the 
10  load  factor  variables  are  bounded  above  and  below. 

TABLE  2. 

Typical  Results  for  a  Ten  Step  Problem 


Program  Outer  Iterations  Total  Function  Execution  Time 

name  performed  Evaluations  (double  precision) 


GRG 

61 

1555 

102  sec 

LPNLP 

262 

6026 

238  sec 

VF01AD 

5 

5481 

289  sec 

VF02AD 

11 

231 

302  sec 

It  should  be 

noted  that 

the  GRG  program  consistently  performed  better 

than  the  others. 

The  LPNLP 

program  was  somewhat  sensitive 

in  that  it 

sometimes  did  converge  taking  approximately  1J  the  time  of  GRG,  while  at 


i 
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other  times  very  long  runs  resulted  with  no  convergence  obtained. 

The  programs  VF01AD  and  VF02AD  were  more  reliable  than  IPNLP,  but 
the  solution  time  and  reliability  of  these  did  not  surpass  GRG. 

The  weighting  factor  in  the  cost  functional  should  make  a  difference 
in  the  optimal  trajectory  calculated.  Weighting  deviations  from  the 
desired  yd(k)  heavily  should  force  the  trajectory  to  stay  close  to  yd, 
at  a  cost  of  increased  altitude.  In  Figure  14  we  illustrate  the  effect 
of  weighting  on  the  solution  of  the  problem. 


cross  track 

yd  =  1000 

(a)  -  weight  =  1.0 

(b)  -  weight  =  5.0 


distance  along  track 

( a )  -  weight  =  1.0 

(b)  -  weight  =  5.0 

(c)  -  terrain  height  along  the 

path  yd 


Fig.  14b 


Fig.  14a 


altitude(ft)  altitude(ft) 
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ALTITUDE  OF  THE  FLIGHT  PATH  (a)  AND  TERRAIN  HEIGHT  (b) 


FOR  THE  CASE  OF  WEIGHT  =1.0 
Fia.  14c 


distance 
between 
flight 
&  terrain 


FOR  THE  CASE  OF  WEIGHT  =  5.0 


Fig  14d 


Notice  that  a  lower  altitude  if  possible  when  the  flight  path  Is 
allowed  to  go  around  the  obstacle. 
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After  working  with  the  special  case  of  stylized  continuous  terrain, 
as  in  the  previous  discussion,  we  incorporated  the  terrain  generation 
model  that  was  described  in  Section  III.  The  computational  results 
were  similar  to  what  had  already  been  observed.  The  6R6  method  was 
the  most  successful  of  the  four.  LPNLP  was  again  erratic  in  its 
performance.  In  general  the  solution  times  were  slightly  less  than  was 
shown  in  Table  2  for  cases  where  the  terrain  was  fairly  smooth,  while 
comparable  to  the  values  in  Table  2  for  rougher  terrain.  The  VF02AD 
and  GRG  programs  were  more  reliable  at  finding  a  better  (local)  minimum 
than  were  either  of  the  augmented  Lagrangian  methods.  There  were 
convergence  problems  in  several  cases  of  the  LPNLP  program. 

Another  factor  that  seemed  to  have  minor  effect  was  the  choice  of 
approximation  method  for  non-grid  terrain  points.  No  consistent  trends 
could  be  seen  in  comparing  the  computational  performance  of  the  bilinear, 
bicubic  Hermite,  and  the  convolutional  cubic  methods.  We  suspect  that 
this  is  so  partly  because  the  calculation  of  the  terrain  height  and 
terrain  slopes  were  a  minor  part  of  the  total  function  evaluation, 
providing  only  ten  of  the  134  total  constraints.  The  discontinuity 
of  the  derivative  of  the  bilinear  method  would  also  be  unlikely  to  arise 
in  practice. 

Problem  Simplifications 

It  is  apparent  from  the  calculation  times  of  Table  2  that  none  of 
the  programs  is  capable  of  finding  an  optimal  ten  segment  path  in  any 
time  approaching  real  time.  For  example,  it  would  be  desired  that  a 
ten  second  segment  of  flight  path  be  recalculated  each  second.  We 
attempted  to  modify  the  statement  of  the  problem  in  order  to 
make  it  somewhat  more  tractable,  while  still  maintaining  its  basic 


structure. 


One  such  simplification  was  to  allow  a  non-uniform  time  step  in 
the  control  sequence.  As  the  problem  has  been  stated  the  total  time  t^ 
was  broken  into  constant  control  intervals  h  seconds  in  length.  Since 
only  the  first  part  (an  amount  tu  seconds)  of  the  total  tf  seconds  of 
flight  path  to  be  computed  will  be  actually  flown  before  a  new  path 
calculation  is  performed,  there  may  be  no  need  to  treat  all  time  steps  as 
of  equal  length.  For  example,  a  20  second  period  could  be  done  in  steps 
of  lengths 

1,  1,  1.5,  3.0,  5.0,  8.5  . 

The  twenty  second  period  would  therefore  be  covered  in  only  six  steps 
instead  of  twenty.  As  long  as  the  constraints  for  all  portions  of  the 
trajectory  are  satisfied,  a  major  saving  in  computation  could  be  realized. 

If  the  initial  portions  of  the  trajectory  of  the  problem  for  six  non-uniform 
steps  turned  out  to  be  similar  to  the  uniform  step,  twenty  segment 
solution  there  would  be  reason  to  use  the  nonuniform  step  method.  The 
optimal  path  generated  by  the  nonuniform  step  method  would  be  less  good 
than  the  full  method,  but  this  trajectory  would  be  safe  and  within  the 
capability  of  the  aircraft.  This  is  because  the  demand  for  a  longer 
period  between  control  adjustment  causes  the  algorithm  to  become  conserv¬ 
ative. 

An  example  of  the  above  procedure  was  run.  A  six  step  problem 
was  made  to  cover  10  seconds  by  choosing  the  intervals  as 
1,  1.1,  1.3,  1.7,  2.3,  2.6  . 

The  time  for  solution  of  the  ten  step  problem  was  314  seconds  using  the 
VF01AD  procedure,  while  the  same  method  applied  to  the  modified  problem 
reached  a  solution  in  89.  seconds.  The  trajectories  obtained  were 
comparable,  with  the  non-uniform  (i.  e.  larger  average  step)  case  taking 
a  wider  cross-track  swing  due  to  the  inability  of  the  controls  of  that 
"airframe"  to  respond  as  quickly. 
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Weighting  early  values  in  the  cost  functional  more  than  later  ones 
would  presumably  force  the  algorithm  to  make  the  initial  part  of  the  path 
approach  its  optimum  in  fewer  iterations  than  would  otherwise  be  the 
case.  Thus,  it  would  be  conceivable  that  the  optimization  could  be  stopped 
sooner.  Provided  that  all  constraints  are  satisfied,  this  suboptimal 
trajectory  might  be  satisfactory.  The  weighting  scheme  discussed  above 
would  change  the  cost  functional  of  equation  (6)  to  the  form 

£  QiPz(i)2  +  W.  [  py(i )  -  yd(i)  ]2 

where  Qi  would  be  larger  for  the  first  several  values  than  for  subsequent 
ones.  We  implemented  such  a  procedure,  and  also  allowed  for  zero  weights 
in  the  cost  functional  for  the  later  segments.  Several  simulations  of 
a  procedure  such  as  this  indicated  that  computation  time  was  not  decreased, 
and  that  the  trajectories  obtained  were  significantly  different  from  those 
of  the  original  implementation. 

Finally,  it  would  be  true  that  in  practice  the  solution  to  the 
previous  N  segment  problem  would  be  provided  as  an  initial  guess  for  the 
next  problem.  Our  results  showed  that  a  considerable  savings  in  computa¬ 
tion  time  may  be  realized,  but  even  with  this,  solution  times  are 
considered  to  be  too  long  for  the  proposed  application.  In  a  typical 
case  a  savings  of  approximately  50%  in  overall  computation  time  would 
occur  when  a  good  initial  guess  was  made  available. 
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VII.  CONCLUSION 

We  set  out  at  the  start  of  this  short  project  to  attempt  to  make 
a  contribution  to  the  general  field  of  nonlinear  programming  algorithms, 
as  well  as  to  investigate  the  properties  of  a  particular  problem. 
Unfortunately,  several  ideas  did  not  come  to  fruition  as  we  had  hoped. 

For  example,  during  the  course  of  the  project  the  author  tried  to 
understand  and  demonstrate  the  practical  consequences  of  the  theoretical 
results  on  parallel  dynamic  programming  as  stated  in  the  paper  by 
Bertsekas  (ref.  47).  Little  progress  was  made  in  the  implementation 
of  this  material  and  there  are  no  results  to  report. 

In  a  similar  vein,  the  use  of  parallel  processing  in  problems  of 
the  sort  arising  in  nonlinear  programming  has  yet  to  be  accomplished. 

The  communications  aspect  and  coordination  of  search  data  has  been  a 
stumbling  block. 

What  we  have  accomplished  is  to  add  to  the  body  of  user  information 
with  regard  to  the  current  state  of  research  in  nonlinear  programming 
algorithms.  We  selected  and  tested  a  set  of  programs  that  illustrate 
the  major  trends  in  the  field,  and  have  shown  that  in  fact  these  methods 
are  successful . 

We  pointed  out,  however,  that  the  current  methods  do  not  seem 
capable  of  meeting  the  computational  demands  of  the  terrain  following 
and  terrain  avoidance  path  generation  problem.  Several  possible 
simplifications  and  modifications  of  the  problem  statement  could  be 
Introduced  in  order  to  speed  solution. 

The  author  is  grateful  to  the  Air  Force  Office  of  Scientific 
Research  for  providing  mini  grant  support  for  the  accomplishment  of 
the  project.  Having  been  involved  in  this  study,  the  author  intends  to 
pursue  research  into  some  of  the  unsolved  areas  that  he  was  exposed  to. 
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